{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian Process Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Distribution\n",
    "![normal dist](./normal.png)\n",
    "- The multivariate **Gaussian** (or **Normal**) distribution has a joint probability density given by\n",
    "$$p(\\mathbf{x} | \\mathbf{m}, \\Sigma)=(2 \\pi)^{-D / 2}|\\Sigma|^{-1 / 2} \\exp \\left(-\\frac{1}{2}(\\mathbf{x}-\\mathbf{m})^{\\top} \\Sigma^{-1}(\\mathbf{x}-\\mathbf{m})\\right)$$\n",
    "where $\\mathbf{m}$ is the mean vector (of length D) and $\\Sigma$ is the (symmetric, positive definite) covariance matrix (of size $D \\times D )$. As a shorthand we write $\\mathbf{x} \\sim \\mathcal{N}(\\mathbf{m}, \\Sigma)$\n",
    "- Let $x$ and $y$ be jointly Gaussian random vectors\n",
    "$$\n",
    "\\left[\\begin{array}{l}{\\mathbf{x}} \\\\ {\\mathbf{y}}\\end{array}\\right] \\sim \\mathcal{N}\\left(\\left[\\begin{array}{c}{\\boldsymbol{\\mu}_{x}} \\\\ {\\boldsymbol{\\mu}_{y}}\\end{array}\\right],\\left[\\begin{array}{cc}{A} & {C} \\\\ {C^{\\top}} & {B}\\end{array}\\right]\\right)=\\mathcal{N}\\left(\\left[\\begin{array}{c}{\\boldsymbol{\\mu}_{x}} \\\\ {\\boldsymbol{\\mu}_{y}}\\end{array}\\right],\\left[\\begin{array}{cc}{\\tilde{A}} & {\\tilde{C}} \\\\ {\\tilde{C}^{\\top}} & {\\tilde{B}}\\end{array}\\right]^{-1}\\right)\n",
    "$$\n",
    "then the marginal distribution of $\\mathbf{x}$ and the conditional distribution of $\\mathbf{x}$ given $\\mathbf{y}$ are\n",
    "$$\n",
    "\\mathbf{x} \\sim \\mathcal{N}\\left(\\boldsymbol{\\mu}_{x}, A\\right), \\text { and } \\mathbf{x} | \\mathbf{y} \\sim \\mathcal{N}\\left(\\boldsymbol{\\mu}_{x}+C B^{-1}\\left(\\mathbf{y}-\\boldsymbol{\\mu}_{y}\\right), A-C B^{-1} C^{\\top}\\right)\n",
    "$$\n",
    "or $$\n",
    "\\mathbf{x} | \\mathbf{y} \\sim \\mathcal{N}\\left(\\boldsymbol{\\mu}_{x}-\\tilde{A}^{-1} \\tilde{C}\\left(\\mathbf{y}-\\boldsymbol{\\mu}_{y}\\right), \\tilde{A}^{-1}\\right)\n",
    "$$\n",
    "- The product of two Gaussian random variables is **NOT** Gaussian, but the product of two Gaussian densities gives another (un-normalized) Gaussian\n",
    "$$\n",
    "\\mathcal{N}(\\mathbf{x} | \\mathbf{a}, A) \\mathcal{N}(\\mathbf{x} | \\mathbf{b}, B)=Z^{-1} \\mathcal{N}(\\mathbf{x} | \\mathbf{c}, C)\n",
    "$$\n",
    "$$\n",
    "\\mathbf{c}=C\\left(A^{-1} \\mathbf{a}+B^{-1} \\mathbf{b}\\right) \\text { and } C=\\left(A^{-1}+B^{-1}\\right)^{-1}\n",
    "$$\n",
    "$$\n",
    "Z^{-1}=(2 \\pi)^{-D / 2}|A+B|^{-1 / 2} \\exp \\left(-\\frac{1}{2}(\\mathbf{a}-\\mathbf{b})^{\\top}(A+B)^{-1}(\\mathbf{a}-\\mathbf{b})\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Gaussian random variable sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# standard Gaussian\n",
    "\n",
    "x = np.random.randn()\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# multivariate Gaussian with covariance matrix C\n",
    "\n",
    "m = np.array([0, 0])\n",
    "C = np.eye(2)\n",
    "D = 1000\n",
    "x = np.random.multivariate_normal(mean=m, cov=C, size=D)\n",
    "plt.plot(x[:, 0], x[:, 1], '.')\n",
    "plt.axis('equal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "x = stats.multivariate_normal(mean=m, cov=C)\n",
    "x.rvs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian Process\n",
    "\n",
    "- A GP is determined by a mean function and covariance function. \n",
    "$$\n",
    "f(\\mathbf{x}) \\sim \\mathcal{G P}\\left(m(\\mathbf{x}), k\\left(\\mathbf{x}, \\mathbf{x}^{\\prime}\\right)\\right)\n",
    "$$\n",
    "$$\n",
    "\\begin{aligned} m(\\mathbf{x}) &=\\mathbb{E}[f(\\mathbf{x})] \\\\ k\\left(\\mathbf{x}, \\mathbf{x}^{\\prime}\\right) &=\\mathbb{E}\\left[(f(\\mathbf{x})-m(\\mathbf{x}))\\left(f\\left(\\mathbf{x}^{\\prime}\\right)-m\\left(\\mathbf{x}^{\\prime}\\right)\\right)\\right] \\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Covariance function or **kernel**\n",
    "![cov](./cov_fun.png)\n",
    "where $d = x-x'$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Generate gaussian process"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cov_se(x, x_, l):\n",
    "    '''\n",
    "    Sequare Exponential kernel.\n",
    "    '''\n",
    "    return np.exp(-np.sum((x - x_)**2) / 2 / l**2)\n",
    "def gen_cov(x, cov_fun):\n",
    "    '''\n",
    "    Generate covariance matrix.\n",
    "    '''\n",
    "    D = len(x)\n",
    "    K = np.empty((D, D))\n",
    "    for i in range(D):\n",
    "        for j in range(D):\n",
    "            K[i,j] = cov_fun(x[i], x[j])\n",
    "    return K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "\n",
    "D = 50\n",
    "n_sample = 5\n",
    "x = np.linspace(0, 1, D)\n",
    "cov_fun = partial(cov_se, l=0.1)\n",
    "K = gen_cov(x, cov_fun)\n",
    "\n",
    "# generate samples of a GP with SE covariance function\n",
    "for i in range(n_sample):    \n",
    "    f = np.random.multivariate_normal(mean=np.zeros(D), cov=K)\n",
    "    plt.plot(x, f)\n",
    "plt.title('GP samples')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)');\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Exercise: use the **periodic** covariance function to generate samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write code here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
